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Abstract 


A grid generation and flow solution algorithm for the Eider equations 
on unstructured grids is presented. The grid generation scheme , which 
utilizes Delaunay triangulation , generates the field points for the mesh 
based on cell aspect ratios and allows clustering of grid points near solid 
surfaces. The flow solution method is an implicit algorithm in which the 
linear set of equations arising at each time step is solved using a Gauss- 
Seidel procedure that is completely vectorizable. In addition , a study 
is conducted to examine the number of subiterations required for good 
convergence of the overall algorithm . Grid generation results are shown 
in two dimensions for an NACA 0012 airfoil as well as a two-element 
configuration. Flow solution results are shown for a two-dimensional 
flow over the NACA 0012 airfoil and for a two-element configuration 
in which the solution has been obtained through an adaptation procedure 
and compared with an exact solution. Preliminary three-dimensional 
results also are shown in which the subsonic flow over a business jet is 
computed . 


Introduction 

The use of unstructured grids for a solution of the 
Euler equations offers several advantages over the use 
of structured grids. These advantages include the 
ease with which adaptive methodology can be incor- 
porated into the flow solvers and the relatively short 
time to generate grids about complex configurations. 
Although the overall time to generate grids about 
complex configurations is much shorter for unstruc- 
tured grids compared with that of block-structured 
grids, the computer time required for the unstruc- 
tured flow solvers historically has been much longer 
than that of structured grids. Although unstruc- 
tured flow solvers will continue to require longer com- 
puter times than those of structured grids because 
of indirect addressing, recent advances (refs. 1 and 
2) make three-dimensional computations on unstruc- 
tured grids more competitive with those of structured 
grids. 

As mentioned, the success of unstructured grids 
mainly is due to the relative ease at which grids 
can be obtained over complex configurations. Two 
dominant methods of generating unstructured grids 
currently exist. The first of these techniques is the 
advancing-front method in which the cells that make 
up the interior of the mesh are computed by marching 
away from the domain boundaries (refs. 3 and 4). 
This method has been used with success to generate 
grids about many complex configurations (ref. 5). 
Further details of this technique can be found in 
references 3 to 6 and the references contained therein. 

The other method commonly used for genera- 
tion of unstructured grids is Delaunay triangulation 


(refs. 7 and 8), which is emphasized in the current 
study. This approach triangulates a given set of 
points in a unique way so that the minimum angle 
of each triangle in the mesh is maximized. The ad- 
vantage of this technique is that the resulting meshes 
are optimal for the given point distribution because 
they do not usually contain many extremely skewed 
cells. 

The field points for generating grids using the 
Delaunay triangulation approach usually arc speci- 
fied a priori by generating points about individual 
components with structured grids (ref. 9), by sub- 
dividing existing quadrilateral cells using a quadtree 
encoding method (ref. 6), or by embedding the ge- 
ometry into a Cartesian grid (ref. 10). A novel ap- 
proach to the generation of field points is given by 
Holmes and Snyder (ref. 11); in this approach, the 
field points are generated as the triangulation pro- 
ceeds based on the aspect ratio and cell area of cur- 
rent triangles. This technique generates grids that 
are not highly skewed because new points are intro- 
duced to continually reduce the cell aspect ratios. 
Unfortunately, grids generated in this manner gener- 
ally are too coarse to be used for obtaining accurate 
flow- field solutions without adaptation. 

In the present study, an approach similar to that 
of Holmes and Snyder is used, and an extension is 
incorporated which automatically adds new nodes 
to cluster points in the regions of interest. Using 
the new generator, grids that are suitable for com- 
putations are efficiently generated around complex, 
multibody configurations. 


Many advances also have been made in flow 
solvers for obtaining flow-field solutions on unstruc- 
tured grids. Impressive results, in which solutions are 
found for a wing configuration using a node-based, 
central-differencing scheme with multigrid to achieve 
rapid convergence, have been obtained by Mavriplis 
(ref. 2). In this reference, solutions on a three- 
dimensional grid consisting of more than 2 million 
cells are obtained in approximately 1 hour. 

For upwind solvers, Frink et al. (ref. 1) have gen- 
erated results for many steady-state applications us- 
ing a cell-centered, multistage time-stepping scheme 
and Roc’s approximate Riemann solver (ref. 12). For 
unsteady applications, Batina (ref. 13) has developed 
both explicit and implicit algorithms for obtaining 
aeroelastic applications, and Rausch et al, (ref. 14) 
have coupled some of these methods with adaptive 
mesh refinement . 

In this study, an implicit algorithm for solving the 
Euler equations is described. This method, along 
with the work in references 13 and 15, is based on 
the backward-Euler time-differencing scheme, but it 
is formulated in a manner that permits full vectoriza- 
tion. In addition, the number of subiterations nec- 
essary to sufficiently solve the linear problem and to 
obtain the best convergence rate is examined. Re- 
sults are shown for both two- and three-dimensional 
calculations. 

The author acknowledges Daryl Bonhaus for 
generating the grid around the business jet. 

Symbols 

A matrix 

A area of cell 

a speed of sound 

body conditions on body 

CFL Cour ant- Friedrichs- Lewy number 

C p pressure coefficient 

c chord length 

D diagonal components of A 

d distance to nearest surface node 

E total energy per unit volume 

F fluxes of mass, momentum, and 

energy 

F fluxes normal to cell face 

split fluxes 
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function used for clustering grid 
points 

identity matrix 
length of cell face i 
components of A below diagonal 
Mach number normal to cell face 
frec-stream Mach number 
components of A above diagonal 
total number of cells 
unit normal 

number of edges meeting at node 
x and y components of unit normal 
all off-diagonal components of A 
pressure 

conserved state vector, 

Q = [p pu pv E] t 

primitive state vector, 
q = [p u v p] T 

residual for cell 

Riemann invariants 

vector from center of cell to center 
of edge 

reference condition 

entropy 

time 

velocity normal to cell face 

Cartesian velocities in x and y 
directions 

Cartesian coordinates 
angle of attack 

parameter used for grid clustering 

ratio of specific heats, taken as 1.4 

percent of spanwise location on 
wing 

density 

standard deviation of $ 
function used for grid clustering 
average value of 
boundary of cell 
relaxation factor 
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Two-Dimensional Grid Generation 

Delaunay Triangulation 

The foundation of the proposed grid generation 
procedure is the Delaunay triangulation method de- 
scribed in detail in reference 7. In this technique, a 
set of points is triangulated by inserting each point, 
one at a time, into a current triangulation so that 
no vertex from one triangle lies within the circumcir- 
cle of any other triangle. The procedure is initiated 
by first identifying all the cells that have a circum- 
circle enclosing the point to be inserted. An exam- 
ple is shown in figure 1; in this figure, the point to 
be inserted lies within the circumcircle of two trian- 
gles. The Delaunay cavity, shown in figure 2, then 
is formed from the union of all the triangles identi- 
fied previously. At this stage, a new triangulation 
is made by simply connecting the new point to each 
of the nodes lying on the boundary of the Delaunay 
cavity, as depicted in figure 3. 



Figure 1. Identifying cells broken by introducing new point. 

To generate grids about arbitrary two- 
dimensional configurations, an initial triangulation 
consisting of a square divided into two triangles is 
formed first. This square has four corner points lo- 
cated a sufficient distance from all solid surfaces. The 
points that define the solid surfaces then are inserted 
using Bowyer’s algorithm (ref. 7), followed by a pre- 
determined number of far-field points that are lo- 
cated in a circular pattern which is a specified radius 
from the center of the bodies. The cells that make up 
the interior of the body then are identified according 



Figure 2. Delaunay cavity. 



Figure 3. Reconnection of grid after inserting new point. 

to whether the center of each cell is located inside or 
outside one of the bodies. 

After this initial phase of the process, a loop is 
conducted over all the cells, and a new point is imme- 
diately introduced at the center of the circumcircle of 
any triangle that has an aspect ratio (defined as the 
ratio of the circumcircle radius to twice the in-circle 
radius) exceeding a predetermined tolerance of ap- 
proximately 1.5. The surface integrity is maintained 
by rejecting any point that would result in breaking 
the cells that make up the airfoil interior (ref. 8). 
Note that when a cell aspect ratio is larger than the 
tolerance, the new point is immediately added into 
the existing triangulation. This addition prevents 
duplicate points from being added when two triangles 
have points that define the same circumcircle. 
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Immediately adding points in this manner elim- 
inates the searching that is otherwise necessary to 
identify the first triangle, which is broken by the ad- 
dition of the current point. The elimination is possi- 
ble because the cell that has an aspect ratio greater 
than the tolerance will also correspond to one of the 
triangles which forms the Delaunay cavity. Because 
no searching is required, the computer time necessary 
to generate the field points in this manner is small. 
This addition of new points is similar to that of the 
method used in reference 1 1 in which new points are 
introduced based on both the cell area and the aspect 
ratio. 

An example of this process is shown in figures 4 
through 9 for a sample grid around an NACA 0012 
airfoil. The airfoil surface is defined with 200 points 
along the surface and 32 points placed around the 
outer boundary. Note that the outer boundary is 
placed close to the airfoil for illustrative purposes, 
thus allowing the entire grid to be seen. Figure 4 
shows the initial triangulation in which only the 
surface points and the outer boundary points have 
been included; this triangulation has cells with aspect 
ratios as high as 160. 



Figure 4. Initial sample grid around NACA 0012. 

New points now are introduced at the center of 
the circumcircle of any cell that has an aspect ratio 
exceeding 1.5. Figures 5 to 8 show a few of the 
intermediate triangulations after inserting the first, 
second, third, and fourth points, respectively. 



Figure 5. Sample grid around NACA 0012 after inserting one 
point. 



Figure 6. Sample grid around NACA 0012 after inserting two 

points. 

The final grid obtained by adding field points in 
this manner is shown in figure 9. This grid, which 
consists of 1328 nodes, 3748 faces, and 2420 cells, 
has a maximum aspect ratio of 1.495. Although all 
the resulting cells are nearly equilateral, the grids 
generated with this technique are coarse a short 
distance from the airfoil and arc not sufficient for 
accurate computations. Therefore, increasing the 
grid density in the vicinity of the airfoil is necessary. 
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Figure 7. Sample grid around NACA 0012 after inserting 
three points. 



Figure 8. Sample grid around NACA 0012 after inserting four 
points. 


Extensions for Clustering Mesh Points 

To add new points in the vicinity of the airfoil, 
a value is first assigned to each existing cell; this 
value is the product of the cell area and a weighting 
function that decreases as the distance from the cell 
center to a solid surface increases: 

4>{A, d) = A x f(d) (1) 

In this equation, d is the distance from the cell center 
to the nearest node that lies on a solid boundary. 


Figure 9. Final sample grid around NACA 0012 with all 
aspect ratios <1.5. 


This variable will be used to add subsequent points 
in cells in which the deviation of from the average is 
larger than the standard deviation. For this reason, 
the average and standard deviations of this variable 
are first computed: 


~4> = 


i N 
— 

N 4-f 

7 = 1 


( 2 ) 


<7 = 


N 


N 


E - fa ) 2 

i=l 

N 


( 3 ) 


A list of new points that will be inserted into the ex- 
isting grid then is constructed from the cell centers 
of all triangles in which the local value of <p(A,d) 
exceeds that of the average plus the standard devi- 
ation, i.c., whenever (j>i > <j> + cr. This list of new 
points then is introduced as before, using Bowyer’s 
algorithm. By adding new points in this manner, the 
function (f> tends to be evenly distributed over the 
grid, and new points are added first at larger cells 
near the body. Few, if any, new points are introduced 
far from the solid surfaces. 

The weighting function used in the current study 
is given by 

(4> 

In this equation, is a distance that is measured 
from the airfoil surface; clustering will occur 




predominantly in regions where the distance to the 
airfoil surface is less than d$. A plot of this func- 
tion is shown in figure 10 for several values of /? 
and do = 0.5. As seen, the transition of this func- 
tion at d — 0.5 steepens as /? increases, and the value 
decreases as the distance from the airfoil increases. 
Thus, the transition between clustered and nonclus- 
tered regions can be made smoothly, and the dis- 
tance away from the airfoil in which clustering occurs 
is also controlled. Note that because this procedure 
only adds one point at the center of each triangle, the 
amount of clustering for the final grid (i.c., how many 
new points are introduced) is increased by repeat- 
ing this procedure several times. In practice, three 
or four repetitions in which (3 is gradually increased 
lead to grids with good clustering near the surface 
of the airfoils, and a reasonably smooth transition 
region between the clustered and nonclustered areas 
is obtained. Further enhancements to this procedure 
may be achieved by varying the weighting function. 


16 r 3=5 



Figure 10. Weighting function for several values of /? and 
= 0.5. 

The final step is to smooth the grid with a simple 
Laplacian-type procedure as given in reference 16. 
This process is achieved by repositioning the mesh 
points according to 

*" +1 

k=l 

n } ( 5 ) 

y? +l = y? + % £ to-vi) 

k=l > 

where u) is a relaxation factor and the sum is obtained 
over all edges meeting at node i. For the current 


study, a relaxation factor of 0.2 is typically used, and 
100 to 200 iterations of smoothing are performed. 

The final sample grid for the NACA 0012 airfoil 
is shown in figure 11. This grid, which demonstrates 
the success of the clustering procedure, is a clear 
improvement to the grid previously shown in figure 9. 



Figure 11. Final sample grid around NACA 0012. 


Although the techniques just outlined currently 
have not been implemented in three dimensions, no 
readily apparent implementation obstacle exists. 

Euler Solver 

The Euler flow solver is an implicit, cell-centered, 
upwind-differencing code in which the fluxes on the 
cell faces are obtained using the Van Leer flux- 
vector splitting technique (ref. 17). The solution 
at each time step is updated using an implicit algo- 
rithm that uses the linearized, backward- Euler, time- 
differencing scheme. At each time step, the linear 
system of equations is solved with a subiterative pro- 
cedure in which the mesh cells are divided into groups 
(or colors) so that no two cells in a given group share 
a common edge. For each subiteration, the solution 
is obtained by solving all the unknowns in a given 
color before proceeding to the next color. Because 
the solution of the unknowns in each gToup depends 
on those from previous groups, a Gauss-Seidel-type 
procedure that is completely vectorizable is obtained. 

Governing Equations 

The governing equations are the time-dependent 
Euler equations, which express the conservation of 
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mass, momentum, and energy for an inviscid gas. 
The equations are given by 


dQ 1 


j> FhdQ = 

n 


0 


are given as 
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F n 


P 

pu 

pv 

E 


pU 

pUu + h x p 

pUv + TLyP 

( E + p)U 


(8) 


and U is the velocity in the direction of the outward 
pointing unit normal to a cell face 


U = h x u + hyV 


(9) 


The equations are closed with the equation of state 
for a perfect gas 


P = 


(7-1 )[E-p (u 2 + v 2 ) /2 


( 10 ) 


the Jacobian matrix of F has negative eigenvalues. 
The split fluxes are given by 


(6) 

F ± = 

n 

where 

(7) 

and 




/ mass 

/mass {fax (~U ± 2a) h\ + a} 
/mass { [% (-U ± 2a) + t'} 

/*± 

/energy 


/mass = ( M n ± 1) / 4 


} (13) 


(14) 


r± - f± 

/energy — /mass 


(1 -7)J7 2 ±2(7- \)Ua + 2a 2 u 2 + v 2 

7 2 _ i 2 


The steady-state residual, given by 
R = — j) F • n dQ 


(15) 

(16) 


is calculated using a trapezoidal integration by sum- 
ming the fluxes over each of the faces that make up 
the control volume. For example, the residual in a 
triangular cell is calculated as 


f - i= 3 

R = - j> F.ndn = -^[F+ (Q-) 


+ F- (Qt)]i f (17) 


Flux- Vector and Residual Calculation 

For the computations shown in this report, the 
flux vectors in equation (8) are upwind differenced 
using the flux- vector splitting technique of Van Leer 
(ref. 17). These flux vectors are given in terms of 
the Mach number normal to the cell face, defined as 
M n = Ufa. For supersonic flow in the direction of a 
face normal ( M n > 1) 

F + = (f • + = F F _ = (f • n) -0(11) 

whereas for supersonic flow in the opposite direction 
of the face normal (M n < —1) 

F- = (F n) =F F+ = ^nj + =0 (12) 

For subsonic flow ( | M n | < 1), the fluxes are split 
into two contributions, F + and F _ , such that the 
Jacobian matrix of F + has positive eigenvalues and 


Here F ± (Q =I= ) represents the split fluxes on the cell 
faces formed from an upwind interpolation of the 
data to each face. For first-order accurate differenc- 
ing, the data on the face arc obtained from the data 
in the cells that lie on each side of the cell face. For 
higher order differencing, the primitive variables arc 
extrapolated to the cell faces using a Taylor scries ex- 
pansion about the center of the cell so that the data 
on the face are given by 

Qface = ^center d" VQ ' r (^®) 

where r is the vector extending from the center of 
the cell to the center of the cell face. 

For evaluating the gradient the data first 

are interpolated to the nodes using inverse distance 
weighting, and the gradient then is evaluated using 
Green’s theorem. This interpolation method is fur- 
ther discussed in reference 18. Note that obtaining 
the data at the nodes also has been accomplished 
using a linear least-squares fit of the data in the sur- 
rounding cells with no apparent differences observed 
in the solutions obtained with either method. 
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Boundary Conditions 

The boundary conditions on the body are set ac- 
cording to characteristic-type boundary conditions 
similar to those in reference 19. The density, pres- 
sure, and velocity components on the body are 
determined according to 

Pbody = Pref T Pref G ref T fit fV) (19) 

Pbody Pref T (pbody — Pref) / a ref (20) 

w body = u ref ~ n x (n x u + h y v ) ref (21) 

u body = v ref ~ (% u + n y v) rcf (22) 

from which the energy is set using the equation of 
state given in equation (10). The reference conditions 
for equations (19) through (22) are taken from the 
first cell in the grid interior. 

Because an implicit scheme is used in this study, 
implicit boundary conditions are implemented by 
assuming that 


The Cartesian velocities are determined on the 
outer boundary by decomposing the normal and 
tangential velocity vectors into components that 
yield 


^boundary — u ref + (^boundary — ^ref) 
^boundary ~ ^ref T fly (^boundary — ^ref) 

where the subscript ref represents values obtained 
from one point outside the domain for inflow and 
from one point inside the domain for outflow. 

The entropy is determined by using the value 
from either outside or inside the domain, depending 
on whether the boundary is an inflow or outflow 
boundary. Once the entropy is known, the density on 
the far-field boundary is calculated from the entropy 
and speed of sound as 



Pbonndary — 


( ^boundary \ 
7 ^boundary ) 


l 



(31) 


^ Pbody — ^Pref (23) 

A (/™)budy = A (/™)ref “ ( r 7r A if™) + (/”’)) rcf (24) 

A (P°) body = A (Href ~ (fir A iP u ) + fi/ A ^) rc f (^) 

^^body = ^^ref ( 26 ) 

In this manner, the matrix entries that correspond 
to cells lying adjacent to a solid surface can be easily 
modified to include the boundary influence. 

For the far field, explicit boundary conditions are 
used in which the velocity and speed of sound are 
obtained from two locally one-dimensional Riemann 
invariants given by 

R ± = U± (27) 

7 — 1 

These invariants are considered constant along char- 
acteristics that are defined normal to the outer 
boundary. For subsonic conditions at the bound- 
ary, R~ can be evaluated locally from free-stream 
conditions outside the computational domain, and 
R + is evaluated locally from the interior of the do- 
main. The local normal velocity and speed of sound 
on the boundary are calculated using the Riemann 
invariants as 

^boundary = ^ (^ + + R ) (28) 

^boundary “ ^ — R ) (29) 


The energy then is calculated from the equation of 
state. 

Time Advancement Schemes 

Implicit algorithms . The starting point for 
the time advancement algorithm is the linearized, 
backward- Euler, time-differencing scheme that yields 
a system of linear equations for the solution at each 
step given by 


where 


[A] n {AQ}" = {R}" 

(32) 

[Al n — —l + 

1 J A< + dQ 

(33) 


The solution of equation (32) can, in principle, be 
obtained by a direct inversion of [A] 77 ; this solution 
has the advantage of a resulting scheme that becomes 
a Newton iteration in the limit as the time step 
approaches infinity if the exact linearization of R 77 
is used in forming [A] n . Although this technique 
is quite successful in two dimensions (ref. 20), the 
solution at each time step requires a great deal of 
memory to store the components of [A] n as w r ell 
as extensive computer time to perform the matrix 
inversions. This approach, therefore, is currently 
not very feasible for practical calculations in three 
dimensions. 


Because the number of operations required to 
invert a matrix depends on the matrix bandwidth, 
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first-order accurate approximations on the left-hand 
side of equation (32) are often utilized to reduce 
both required storage and computer time. With this 
simplification, the consistency between the left- and 
right-hand sides of equation (32) requires that first- 
order approximations also be used on the right-hand 
side to achieve quadratic convergence. However, with 
first-order approximations on the left-hand (implicit) 
side and second-order approximations on the right- 
hand side, this scheme remains stable for large time 
steps. First-order differencing of the left-hand side 
with higher order differencing on the right-hand side, 
therefore, is considered in the present study. 

A sample configuration of triangles in which the 
cells are randomly ordered is shown in figure 12. The 
corresponding form of the matrix [A]” is shown in 
figure 13; in this figure a circle represents the nonzero 
entries. 
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Figure 12. Sample cell configuration. 


matrices representing the diagonal, subdiagonal, and 
superdiagonal terms, that is 

[A] n = [D] n + [Mf + [N] n (34) 

The simplest iterative scheme for obtaining a solu- 
tion to the linear system of equations is a Jacobi- 
type method in •which all the off-diagonal terms of 
[A] n {AQ> (i.c., [Mf{AQ} + |N] n {AQ}) are taken 
to the right-hand side of equation (32) and are eval- 
uated using the values of {AQ}' from the previous 
subiteration level i. This scheme can be represented 
as 

[D] n {AQ}* +1 = [{R}" - [M + Nf {AQ}'] 

= [{R} n - [Of {AQ}*] (35) 

The disadvantage of this scheme is that the se- 
quence of Jacobi iterations may converge slowly. To 
accelerate the convergence, a Gauss-Seidel procedure 
may be employed in which values of {AQ} are used 
on the right-hand side of equation (35) as soon as 
they are available. An example of this scheme can 
be written as 

[D] {AQ} J+1 = [{R}” - [M] n {AQ} 7+1 - [N] n {AQ} 7 ] (36) 


o ©# 

© © © 

@ m © 



Figure 13. Form of matrix for cells in figure 12. 


Although the solution of the system of equations 
may be obtained through a direct inversion of [A] , 
as previously mentioned, the need for large memory 
can be circumvented through the use of a variety of 
relaxation schemes. In these schemes, the solution 
of equation (32) is obtained through a sequence 
of iterations in which an approximation of AQ is 
continually refined. 

To facilitate the derivation of these schemes, [A] n 
is first written as a linear combination of three 


where the latest values of {AQ} from the subdiag- 
onal terms are immediately used on the right-hand 
side of the iteration equation. A slight modification 
to this algorithm in which the latest values of {AQ} 
from the superdiagonal terms are used results in a 
similar scheme that is given by 


[D] {AQ} 7+1 = [{R} n - [Mj n {AQ} 7 - [N] 77 {AQ} 7 '* 1 ] (37) 

Another variation of this algorithm can be obtained 
by alternating the use of equation (36) with equa- 
tion (37) so that a symmetric Gauss- Seidel- type 
procedure is obtained. 

Note that the algorithms given by equations (36) 
and (37) can both be implemented by sweeping se- 
quentially through each mesh cell and simply using 
the latest values of {AQ} for all the off-diagonal 
terms that have been taken to the right-hand side. 
This procedure can be represented as 


[D] {AQ}* +1 = 


H-r 

{R} 71 — [O]" { AQ} 


(38) 


i+1 

where Q * is the most recent value of Q; this term 
will be at the subiteration level i + 1 for the cells that 
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have been previously updated and at the level i for 
the cells that remain to be updated. The distinc- 
tion between the algorithms in equations (36) and 
(37) comes from sweeping forward through the cells 
(eq. (36)) or backward through the cells (eq. (37)). 

Two disadvantages of this scheme exist. The 
first disadvantage is that this process is not easily 
vectorized because the solution of each point must 
be obtained before proceeding to the next point. 
The second disadvantage is that although the off- 
diagonal terms may be updated and immediately 
used on the right-hand side, the solution of the 
next unknowm may or may not depend on previously 
determined quantities. For example, as can be seen 
from figure 12, when solving for the second unknown 
using equation (36), the updated value of the solution 
at point 1 is not used; therefore, the solution for point 
2 remains a Jacobi step. 

Note that for structured grids in which the cells 
are ordered in a natural manner (e.g., left to right 
and top to bottom), the latest information will im- 
mediately be used for the calculation of the next un- 
known. This occurs because the ordering of the cells 
produces a banded matrix with terms grouped along 
the diagonal. The fact that the latest-obtained data 
arc not necessarily used for updating information in 
unstructured grids is because of the random ordering 
of the cells. 

An improvement to the scheme just described, 
therefore, can be obtained by simply renumbering the 
cells to group terms along the diagonal of the matrix. 
In this manner, the solution of each point will tend 
to ensure that previously updated information from 
the surrounding cells is used as soon as it is available. 
An example of this is shown in figure 14 where the 
same cells used in figure 12 are simply renumbered 
from bottom to top and left to right. The resulting 
form of the matrix, shown in figure 15, shows that 
the grouping along the diagonal is greatly improved. 
The ordering of the cells in this way should result in 
faster convergence of the linear problem than a ran- 
dom ordering of cells. Although the ordering of the 
cells in this example groups unknowns along the di- 
agonal, other procedures such as the Cuthill-McKee 
method described in reference 21 are more effective 
for general configurations. Again, note that several 
variations of this scheme can be obtained by using 
various combinations of equations (36) and (37). An 
important disadvantage of this scheme, however, is 
that the contribution of the off-diagonal terms to 
the right-hand side of equation (38) still cannot be 
vectorized. 


4 

3 

8 

7 

12 

11 

2 

1 

6 

5 

10 

9 N. 


Figure 14. Sample cells. 



Figure 15. Form of matrix for cells in figure 14. 


The Jacobi, Gauss-Seidel, and symmetric Gauss- 
Seidel schemes just described have all been used 
in practice by various researchers. An example of 
research that used these schemes to solve the Euler 
equations for transonic flow over a circular arc in a 
channel is given in reference 15. In this reference, 
the symmetric Gauss-Seidel scheme was shown to 
exhibit the fastest convergence rate of the three 
schemes. The successful use of a symmetric Gauss- 
Seidel algorithm for transonic flow over airfoils is 
described in reference 22; in this work, grouping the 
unknowns along the diagonal is enhanced by sorting 
them according to the x-coordinate direction. 

Vectorization of Gauss- Seidel, The number- 
ing of cells used in the current study is shown in 
figure 16. The ordering is obtained by grouping cells 
so that no two cells in a given group share a common 
edge. The resulting matrix form for [A] is given in 
figure 17. Note that for this example, only two groups 
are formed; in practice, at most, four groups will 
be formed for two-dimensional calculations and five 
groups are formed for three-dimensional calculations. 
The first group for the present example consists of the 
cells numbered 1 through 6, and the second group 
contains cells numbered 7 through 12. 

The solution scheme, which can be written as 
before using equation (38), is implemented by solving 
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Figure 16. Sample cells. 
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Figure 17. Form of matrix for cells in figure 16. 

for all the unknowns in a group at a time. The 
cells in the first group are solved using a Jacobi- type 
iteration, and the cells in all the subsequent groups 
are obtained by using the most recently updated 
values of {AQ} from the off-diagonal contributions. 
In this way, a Gauss- Seidel- type scheme is obtained 
which is easily implemented and is fully vectorized. 
Note that a symmetric Gauss- Seidel- type procedure 
is not necessary and is not used; stability is achieved 
as long as the matrix [A] maintains block diagonal 
dominance that occurs when first-order differencing 
is used on the implicit side of the equation (ref. 23). 

In these discussions, the exact number of sub- 
iterations required to sufficiently converge the linear 
problem (eq. (32)) has not been specified. The num- 
ber of subiterations used for each global time step 
has been determined through numerical experiments 
that are presented in the results. 


Time Step Calculation 

To enhance the convergence to a steady state, 
local time stepping is used. The time step calculation 
for each cell is given by 


At - CFL . . L (39) 

Vr + v z + a 


where L is a length scale for the cell. This length 
scale is defined as the area of the cell divided by the 
perimeter. 

Results 

Flow-field calculations for several demonstration 
cases are presented. The first case is for an 
NACA 0012 airfoil at a free-stream Mach num- 
ber of 0.8 and an angle of attack of 1.25°. The 
grid has an outer boundary placed approximately 
50 chord lengths away from the body; this grid con- 
sists of 3624 nodes, 7012 cells, and 10 636 faces. A 
near- field view of the grid is shown in figure 18. 



Figure 18. Near-field view of grid around NACA 0012 airfoil. 

The pressure coefficient distribution along the 
airfoil surface is shown in figure 19. As seen, a 
moderately strong shock is captured on the upper 
surface of the airfoil, and a weaker shock is captured 
on the lower surface. Also, note that because a flux 
limiter has not been used for the present calculation, 
an “overshoot” is evident ahead of the upper-surface 
shock. The corresponding Mach number contours for 
this case are shown in figure 20. 

For this calculation, the residual of the continu- 
ity equation has been reduced to “machine zero” in 
approximately 400 global iterations, as seen in fig- 
ure 21. The CFL number began at 50 and was 
linearly ramped to 200 throughout 100 iterations. 
The CFL numbers used for the current calculation 
may not be optimal for the present case, but they 
give reasonably good convergence for a wide range 
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Figure 19. Pressure distribution for NAGA 0012 airfoil. 
M x — 0.8; a — 1.25°. 



Figure 20. Mach number contours for NACA 0012 airfoil. 

Afoc = 0.8; a = 1.25°. 

of test problems and grid densities. The memory re- 
quired corresponds to approximately 180 words per 
cell. For each global iteration, 20 subiterations have 
been used "to solve the linear system each time. This 
process results in a computational rate of approx- 
imately 60 /isec per cell per global time step on a 
CRAY YMP using a single processor. This compu- 
tational rate, however, depends on the number of 


subiterations performed. For the current research, 
this number is based on results of a numerical study 
in which the number of subiterations has been var- 
ied for a wide variety of CFL numbers. A typi- 
cal plot of the computer time required to obtain a 
four-order-of-magnitude reduction in the residual is 
shown in figure 22. This plot clearly indicates that 
15 to 20 subiterations for each global iteration pro- 
duce the fastest convergence rate. A similar study 
has been conducted on other grids and other cases 
with similar results. For this reason, between 15 to 
20 sub- iterations are used for all cases shown in this 
report. Although this number of subiterations has 
proved adequate for the current work, further work 
in optimizing this parameter may prove beneficial. 

4 r 
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Figure 21. Convergence history for NACA 0012 airfoil. 
Afoc = 0.8; a — 1.25°. 


A comparison of convergence rates obtained with 
both implicit and explicit boundary conditions on the 
airfoil surface is shown in figure 23. As seen, the use 
of implicit boundary conditions improves the rate of 
convergence through the first several orders of mag- 
nitude. In addition, the use of explicit boundary con- 
ditions impedes the convergence past approximately 
seven orders of magnitude. This behavior has been 
observed for a variety of other cases that are not pre- 
sented. Note that the use of explicit boundary con- 
ditions seems to lead to a more robust code because 
ramping of the CFL number has not been necessary 
when explicit boundary conditions have been used. 
For most calculations in this study, implicit bound- 
ary conditions have been used after five iterations, 
and the CFL is ramped from 50 to 200. 
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Figure 22. Computer time for four-order-of- magnitude re- 
duction in residual for NACA 0012 airfoil. M 0 0 = 0.8; 
a — 1.25°. 
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Figure 23. Comparison of convergence rates for implicit 
and explicit boundary conditions for NACA 0012 airfoil. 
Moc = 0.8; a = 1.25°. 


The next case presented is for a two-element air- 
foil in which an exact incompressible solution exists 
(ref. 24). The initial grid used for this calculation is 
shown in figure 24; this grid consists of 1556 points, 
2882 cells, and 4439 faces. The main element and flap 
have 100 points each along the surface. Note that for 
this calculation, no clustering of cells has been per- 
formed near the surfaces because the final solution is 
obtained through an adaptation procedure described 
in reference 25. The pressure distribution calculated 


at a free-stream Mach number of 0.2 using this ini 
tial grid is shown in figure 25. As seen, the coarse 
grid yields results that agree poorly with the exact 
solution. 

As previously mentioned, a solution also has been 
obtained by adapting the grid to the solution. Adap- 
tation is achieved by first identifying a list of cells 
that require refinement. New points, which are 
located at the center of each of these cells, then 
are introduced into the existing triangulation using 
Bowyer’s algorithm for the Delaunay triangulation, 
and the solution is interpolated to the new grid for 
use in restarting the solution. Because the present 
flow field does not contain discontinuities, the list 
of new cells is identified by flagging all the cells in 
which the undivided velocity gradient exceeds that 
of the average plus the standard deviation of all the 
cells in the grid (refs. 25 and 26). 



Figure 24. Grid around two-element configuration. 

The final grid, shown in figure 26, consists of 3165 
nodes, 6332 cells, and 9190 faces, with 148 nodes on 
the surface of the main clement and 128 nodes on 
the flap. The pressure distribution obtained on this 
grid is shown in figure 27. The agreement with the 
exact solution is improved over that in figure 25. In 
addition, the calculated lift of 2.026 compares well 
with the exact value of 2.0281 given in reference 24. 

The present algorithm also has been implemented 
in three dimensions with the preliminary results sub- 
sequently shown. The case shown is for a business jet 
at Moc = 0.2 and a = 3°. The grid used for the com- 
putations has been generated using the advancing 
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Figure 25. Pressure distribution over two-element configura- 
tion using initial grid. 



Figure 26. Grid obtained for two-element airfoil after adap- 

tat ion. 

front-type grid generation described in reference 27; 
this grid consists of 27 191 nodes, 144 100 cells, and 
294 109 faces. The surface grid for this computa- 
tion, which consists of 11582 triangles, is shown in 
figure 28. 

This case has been run at a constant CFL number 
of 300 with 15 subiterations; the convergence history 
is shown in figure 29. As seen, the residual is re- 
duced between two and three orders of magnitude 
in 100 global iterations, at which point the conver- 


+ Calculations 
— Exact solution 



x/c 

Figure 27. Solution obtained for two-element airfoil after 
adaptation. 



Figure 28. Surface grid for business jet. 

gence rate degrades. After 200 iterations, a resid- 
ual reduction of slightly greater than three orders of 
magnitude is obtained. This “tailing off” behavior 
has not been observed in two dimensions when im- 
plicit boundary conditions are used, but it may be 
caused by the low free-stream Mach number or the 
close proximity of the outer boundary that extends 
approximately 10 body lengths ahead of and behind 
the airplane (but only about two body lengths above 
and below). Note that the “tailing off 7 of the residual 
may also indicate that the high-frequency errors in 
the scheme have been effectively reduced and that 
the low-frequency errors have begun to dominate. 
The use of a multigrid to rapidly eliminate these 
low-frequency errors should enhance the convergence. 
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Figure 29. Convergence history for Lear jet. Moo = 0.2; 
a = 3.0°. 


A pressure distribution comparison at the 
77 = 0.44 span station is made in figure 30 with the 


method described in reference 1. In figure 30, the re- 
sults referred to as FIJN3D are those of the present 
study and USM3D refers to those obtained using the 
computer code of reference 1. (This computer code 
is an upwind finite- volume code that uses multistage, 
time-stepping and flux-difference splitting.) As seen, 
the comparison between the two codes is reasonably 
close, and the main discrepancies occur at the lead- 
ing edge. These differences are due to slight differ- 
ences in the computation of the boundary fluxes and 
because the computations with USM3D use Roe’s 
flux-diffcrencc splitting (ref. 12 ) instead of flux- vector 
splitting. 

The current implementation of this code in three 
dimensions requires approximately 87 words of main 
memory per cell and about 50 words per face of a 
solid-state device (SSD). The computational rate on 
a CRAY YMP is approximately 140 /x sec per cell 
per iteration, based on 15 subiterations per global 
time step. Note that this timing includes both user 
time and system time; without the use of SSD, 
the computational rate improves to approximately 
92 fi sec per cell per iteration because of a significant 
decrease in system time. 
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Figure 30. Comparison of surface pressure distribution. Moo = 0.2\ a = 3.0°; tj - 0.44. 
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Concluding Remarks 

A two-dimensional grid generation procedure has 
been devised which combines automatic point place- 
ment with Delaunay triangulation to efficiently pro- 
duce good-quality unstructured meshes. This 
method, which uses the Delaunay triangulation al- 
gorithm of Bowyer, is based on the work of Holmes. 
The present algorithm improves the previously cited 
work by allowing the automatic generation of new 
mesh points so that the clustering of points near sur- 
faces is achieved. 

A flow solver that is implicit and can be com- 
pletely vectorized is also developed in both two and 
three dimensions. This scheme is based on backward- 
Euler time differencing; the linear problem arising 
at each step is solved by using several iterations of 
a Gauss-Seidel-type procedure. In this method, the 
unknowns are divided into groups so that no cells in 
a given group share an edge; therefore, all the cells in 
a group are independent of each other so that their 
solutions can be obtained simultaneously. 

The effect of the number of subiterations on the 
convergence rate (based on computer time) is also 
examined. Between 15 to 20 subiterations per global 
time step produce the best results. In addition, 
the use of implicit boundary conditions improves the 
convergence rate of the current algorithm. 

Results are shown for a two-dimensional flow over 
an NACA 0012 airfoil and a two-element airfoil in 
which the solution is obtained with adaptation. For 
the two-element configuration, comparisons arc made 
with the exact solution, and excellent results are 
obtained by adaptation. For three dimensions, the 
calculation of subsonic flow over a business jet is 
demonstrated. 

NASA Langley Research Center 
Hampton, VA 23G65-5225 
February 20, 1992 
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